clear
current_directory = pwd;

num_comps = 20
%parpool('local',num_comps)

rng(17)

avg_literature_ag_shocks = importdata('avg_literature_ag_shocks.mat');
ag_productivity_shocks = avg_literature_ag_shocks';
manu_shocks = importdata('manu_productivity_shocks_2080_2099_CSIRO_climateadaptation.mat')';
manu_shocks = ones(1,length(manu_shocks))+manu_shocks;
serv_shocks = importdata('serv_productivity_shocks_2080_2099_CSIRO_climateadaptation.mat')';
serv_shocks = ones(1,length(serv_shocks))+serv_shocks;

% First, run the baseline simulation
baselinesimulation

% Save all the moments from the baseline simulation
% Import shares and domestic expenditure shares for all countries and bilateral pairs
ag_import_shares_mat_baseline = ag_import_share_mat;
ag_domestic_expenditures_vec_baseline = ag_import_share_mat(logical(eye(n_countries)));
manu_import_shares_mat_baseline = manu_import_share_mat;
manu_domestic_expenditures_vec_baseline = manu_import_share_mat(logical(eye(n_countries)));
ag_trade_mat_baseline = ag_trade_mat;
manu_trade_mat_baseline = manu_trade_mat;

% GDP shares by sector 
ag_GDP_share_baseline = sim_ag_shares_GDP;
manu_GDP_share_baseline = sim_manu_shares_GDP;
serv_GDP_share_baseline = sim_serv_shares_GDP;
ag_exp_share_baseline = exp_share_ag_vec;
manu_exp_share_baseline = exp_share_manu_vec;
serv_exp_share_baseline = exp_share_serv_vec;
ag_net_exports_baseline = ag_GDP_share_baseline-ag_exp_share_baseline;
manu_net_exports_baseline = manu_GDP_share_baseline-manu_exp_share_baseline;

% Aggregate sectoral consumption and price indices by country
C_vec_baseline = C_vec;
P_vec_baseline = P_vec;
P_a_vec_baseline = P_a_vec;
P_m_vec_baseline = P_m_vec;
P_s_vec_baseline = P_s_vec;
C_a_vec_baseline = C_a_vec;
C_m_vec_baseline = C_m_vec;
C_s_vec_baseline = C_s_vec;
income_per_capita_baseline = wage_guesses;
total_income_vec_baseline = total_income_vec;
GDP_baseline = P_a_vec_baseline.*C_a_vec_baseline+P_m_vec_baseline.*C_m_vec_baseline+P_s_vec_baseline.*C_s_vec_baseline;
P_Tornqvist_vec_baseline = P_a_vec_baseline.^(ag_exp_share_baseline).*P_m_vec_baseline.^(manu_exp_share_baseline)...
    .*P_s_vec_baseline.^(serv_exp_share_baseline);
% P_Tornqvist_vec_baseline = P_a_vec_baseline.^(ag_GDP_share_baseline).*P_m_vec_baseline.^(manu_GDP_share_baseline)...
%     .*P_s_vec_baseline.^(serv_GDP_share_baseline);
GDP_Tornqvist_baseline = total_income_vec_baseline./P_Tornqvist_vec_baseline;

%%% NO REALLOCATION COUNTERFACTUAL:
P_vec_baseline_check = [];
for cc=1:n_countries
P_vec_baseline_check = [P_vec_baseline_check, ...
    ((omega_a*P_a_vec(cc)^(1-sigma))^((1-sigma)/epsilon_a)*(exp_share_ag_vec(cc)*wage_guesses(cc)^(1-sigma))^(1-(1-sigma)/epsilon_a)...
    +(omega_m*P_m_vec(cc)^(1-sigma))^((1-sigma)/epsilon_m)*(exp_share_manu_vec(cc)*wage_guesses(cc)^(1-sigma))^(1-(1-sigma)/epsilon_m)...
    +(omega_s*P_s_vec(cc)^(1-sigma))^((1-sigma)/epsilon_s)*(exp_share_serv_vec(cc)*wage_guesses(cc)^(1-sigma))^(1-(1-sigma)/epsilon_s))^(1/(1-sigma))];
end
P_a_naive_counterfactual_vec = P_a_vec_baseline./ag_productivity_shocks;
P_m_naive_counterfactual_vec = P_m_vec_baseline./manu_shocks;
P_s_naive_counterfactual_vec = P_s_vec_baseline./serv_shocks;

P_vec_no_reallocation = [];
for cc=1:n_countries
P_vec_no_reallocation = [P_vec_no_reallocation, ...
    ((omega_a*P_a_naive_counterfactual_vec(cc)^(1-sigma))^((1-sigma)/epsilon_a)*(exp_share_ag_vec(cc)*wage_guesses(cc)^(1-sigma))^(1-(1-sigma)/epsilon_a)...
    +(omega_m*P_m_naive_counterfactual_vec(cc)^(1-sigma))^((1-sigma)/epsilon_m)*(exp_share_manu_vec(cc)*wage_guesses(cc)^(1-sigma))^(1-(1-sigma)/epsilon_m)...
    +(omega_s*P_s_naive_counterfactual_vec(cc)^(1-sigma))^((1-sigma)/epsilon_s)*(exp_share_serv_vec(cc)*wage_guesses(cc)^(1-sigma))^(1-(1-sigma)/epsilon_s))^(1/(1-sigma))];
end
C_vec_no_reallocation = wage_guesses./P_vec_no_reallocation;

Baqaee_Farhi_losses_no_reallocation = ag_productivity_shocks.*ag_GDP_share_baseline+manu_shocks.*manu_GDP_share_baseline+serv_shocks.*serv_GDP_share_baseline;

C_a_vec_naive_counterfactual = (P_a_vec_baseline.*C_a_vec_baseline)./(P_a_naive_counterfactual_vec);
C_m_vec_naive_counterfactual = (P_m_vec_baseline.*C_m_vec_baseline)./(P_m_naive_counterfactual_vec);
C_s_vec_naive_counterfactual = (P_s_vec_baseline.*C_s_vec_baseline)./(P_s_naive_counterfactual_vec);

chain_weighted_GDP_naive_counterfactual = calc_chain_weighted_GDP(P_a_vec_baseline,P_m_vec_baseline,P_s_vec_baseline,...
    C_a_vec_baseline,C_m_vec_baseline,C_s_vec_baseline,P_a_naive_counterfactual_vec,P_m_naive_counterfactual_vec,...
    P_s_naive_counterfactual_vec,C_a_vec_naive_counterfactual,C_m_vec_naive_counterfactual,...
    C_s_vec_naive_counterfactual);

P_Tornqvist_vec_no_reallocation = P_a_naive_counterfactual_vec.^(ag_exp_share_baseline).*P_m_naive_counterfactual_vec.^(manu_exp_share_baseline)...
    .*P_s_naive_counterfactual_vec.^(serv_exp_share_baseline);
% P_Tornqvist_vec_no_reallocation = P_a_naive_counterfactual_vec.^(ag_GDP_share_baseline).*P_m_naive_counterfactual_vec.^(manu_GDP_share_baseline)...
%     .*P_s_naive_counterfactual_vec.^(serv_GDP_share_baseline);
GDP_Tornqvist_no_reallocation = total_income_vec_baseline./P_Tornqvist_vec_no_reallocation;

%%% CLOSED ECONOMY COUNTERFACTUAL:
wage_guesses = income_per_capita_baseline;
P_a_naive_counterfactual_vec = P_a_vec_baseline./ag_productivity_shocks;
P_a_vec = P_a_naive_counterfactual_vec;
P_m_naive_counterfactual_vec = P_m_vec_baseline./manu_shocks;
P_m_vec = P_m_naive_counterfactual_vec;
P_s_naive_counterfactual_vec = P_s_vec_baseline./serv_shocks;
P_s_vec = P_s_naive_counterfactual_vec;


calculate_consumption_shares_parallel
ag_exp_share_closed_economy_counterfactual = exp_share_ag_vec;
manu_exp_share_closed_economy_counterfactual = exp_share_manu_vec;
serv_exp_share_closed_economy_counterfactual = exp_share_serv_vec;
delta_ag_GDP_share_closed_economy = ag_exp_share_closed_economy_counterfactual-ag_exp_share_baseline;
delta_manu_GDP_share_closed_economy = manu_exp_share_closed_economy_counterfactual-manu_exp_share_baseline;
delta_serv_GDP_share_closed_economy = serv_exp_share_closed_economy_counterfactual-serv_exp_share_baseline;
ag_GDP_share_closed_economy_counterfactual = ag_GDP_share_baseline+delta_ag_GDP_share_closed_economy;
manu_GDP_share_closed_economy_counterfactual = manu_GDP_share_baseline+delta_manu_GDP_share_closed_economy;
serv_GDP_share_closed_economy_counterfactual = serv_GDP_share_baseline+delta_serv_GDP_share_closed_economy;
C_vec_closed_economy_counterfactual = C_vec;
P_vec_closed_economy_counterfactual = P_vec;
P_a_vec_closed_economy_counterfactual = P_a_vec;
P_m_vec_closed_economy_counterfactual = P_m_vec;
P_s_vec_closed_economy_counterfactual = P_s_vec;
C_a_vec_closed_economy_counterfactual = C_a_vec;
C_m_vec_closed_economy_counterfactual = C_m_vec;
C_s_vec_closed_economy_counterfactual = C_s_vec;
income_per_capita_closed_economy_counterfactual = wage_guesses;
total_income_vec_closed_economy_counterfactual = total_income_vec;

Baqaee_Farhi_losses_closed_economy_counterfactual = ...
    ag_productivity_shocks.*((ag_GDP_share_baseline+ag_GDP_share_closed_economy_counterfactual)/2);

chain_weighted_GDP_closed_economy_counterfactual = calc_chain_weighted_GDP(P_a_vec_baseline,P_m_vec_baseline,P_s_vec_baseline,...
    C_a_vec_baseline,C_m_vec_baseline,C_s_vec_baseline,P_a_vec_closed_economy_counterfactual,P_m_vec_closed_economy_counterfactual,...
    P_s_vec_closed_economy_counterfactual,C_a_vec_closed_economy_counterfactual,C_m_vec_closed_economy_counterfactual,...
    C_s_vec_closed_economy_counterfactual);

P_Tornqvist_vec_closed_economy = P_a_vec_closed_economy_counterfactual.^((ag_exp_share_baseline+ag_exp_share_closed_economy_counterfactual)/2)...
    .*P_m_vec_closed_economy_counterfactual.^((manu_exp_share_baseline+manu_exp_share_closed_economy_counterfactual)/2)...
    .*P_s_vec_closed_economy_counterfactual.^((serv_exp_share_baseline+serv_exp_share_closed_economy_counterfactual)/2);
% P_Tornqvist_vec_closed_economy = P_a_vec_closed_economy_counterfactual.^((ag_GDP_share_baseline+ag_GDP_share_closed_economy_counterfactual)/2)...
%     .*P_m_vec_closed_economy_counterfactual.^((manu_GDP_share_baseline+manu_GDP_share_closed_economy_counterfactual)/2)...
%     .*P_s_vec_closed_economy_counterfactual.^((serv_GDP_share_baseline+serv_GDP_share_closed_economy_counterfactual)/2);
GDP_Tornqvist_closed_economy = total_income_vec_closed_economy_counterfactual./P_Tornqvist_vec_closed_economy;


%%% OPEN ECONOMY COUNTERFACTUAL
ag_T_vec = ag_T_vec.*ag_productivity_shocks;
manu_T_vec = manu_T_vec.*manu_shocks;
serv_T_vec = serv_T_vec.*serv_shocks;

sim_Frechets
P_a_vec = geomean(ag_consumed_prices_mat);
P_m_vec = geomean(manu_consumed_prices_mat);
P_s_vec = geomean(serv_consumed_prices_mat);

find_equilibrium_parallel

sim_ag_shares_GDP = ag_income_vec./total_income_vec;
sim_manu_shares_GDP = manu_income_vec./total_income_vec;
sim_serv_shares_GDP = serv_income_vec./total_income_vec;

ag_import_shares_mat_open_economy = ag_import_share_mat;
ag_domestic_expenditures_vec_open_economy = ag_import_share_mat(logical(eye(n_countries)));
manu_import_shares_mat_open_economy = manu_import_share_mat;
manu_domestic_expenditures_vec_open_economy = manu_import_share_mat(logical(eye(n_countries)));
ag_trade_mat_open_economy_counterfactual = ag_trade_mat;
manu_trade_mat_open_economy_counterfactual = manu_trade_mat;
ag_GDP_share_open_economy_counterfactual = sim_ag_shares_GDP;
manu_GDP_share_open_economy_counterfactual = sim_manu_shares_GDP;
serv_GDP_share_open_economy_counterfactual = sim_serv_shares_GDP;
ag_exp_share_open_economy_counterfactual = exp_share_ag_vec;
manu_exp_share_open_economy_counterfactual = exp_share_manu_vec;
serv_exp_share_open_economy_counterfactual = exp_share_serv_vec;
ag_net_exports_open_economy_counterfactual = ag_GDP_share_open_economy_counterfactual-ag_exp_share_open_economy_counterfactual;
manu_net_exports_open_economy_counterfactual = manu_GDP_share_open_economy_counterfactual-manu_exp_share_open_economy_counterfactual;
C_vec_open_economy_counterfactual = C_vec;
P_vec_open_economy_counterfactual = P_vec;
P_a_vec_open_economy_counterfactual = P_a_vec;
P_m_vec_open_economy_counterfactual = P_m_vec;
P_s_vec_open_economy_counterfactual = P_s_vec;
C_a_vec_open_economy_counterfactual = C_a_vec;
C_m_vec_open_economy_counterfactual = C_m_vec;
C_s_vec_open_economy_counterfactual = C_s_vec;
income_per_capita_open_economy_counterfactual = wage_guesses;
total_income_vec_open_economy_counterfactual = total_income_vec;

Baqaee_Farhi_losses_open_economy_counterfactual = ...
    ag_productivity_shocks.*((ag_GDP_share_baseline+ag_GDP_share_open_economy_counterfactual)/2);

chain_weighted_GDP_open_economy_counterfactual = calc_chain_weighted_GDP(P_a_vec_baseline,P_m_vec_baseline,P_s_vec_baseline,...
    C_a_vec_baseline,C_m_vec_baseline,C_s_vec_baseline,P_a_vec_open_economy_counterfactual,P_m_vec_open_economy_counterfactual,...
    P_s_vec_open_economy_counterfactual,C_a_vec_open_economy_counterfactual,C_m_vec_open_economy_counterfactual,...
    C_s_vec_open_economy_counterfactual);

P_Tornqvist_vec_open_economy = P_a_vec_open_economy_counterfactual.^((ag_exp_share_baseline+ag_exp_share_open_economy_counterfactual)/2)...
    .*P_m_vec_open_economy_counterfactual.^((manu_exp_share_baseline+manu_exp_share_open_economy_counterfactual)/2)...
    .*P_s_vec_open_economy_counterfactual.^((serv_exp_share_baseline+serv_exp_share_open_economy_counterfactual)/2);
% P_Tornqvist_vec_open_economy = P_a_vec_open_economy_counterfactual.^((ag_GDP_share_baseline+ag_GDP_share_open_economy_counterfactual)/2)...
%     .*P_m_vec_open_economy_counterfactual.^((manu_GDP_share_baseline+manu_GDP_share_open_economy_counterfactual)/2)...
%     .*P_s_vec_open_economy_counterfactual.^((serv_GDP_share_baseline+serv_GDP_share_open_economy_counterfactual)/2);
GDP_Tornqvist_open_economy = total_income_vec_open_economy_counterfactual./P_Tornqvist_vec_open_economy;

% Calculate equivalent variation across scenarios

EV_naive = calc_EV(C_vec_no_reallocation, income_per_capita_baseline, P_a_vec_baseline, P_m_vec_baseline, P_s_vec_baseline, ...
    epsilon_a,epsilon_m,epsilon_s,sigma,omega_a,omega_m,omega_s,num_comps,n_countries);

EV_closed_economy = calc_EV(C_vec_closed_economy_counterfactual, income_per_capita_baseline, P_a_vec_baseline, P_m_vec_baseline, P_s_vec_baseline, ...
    epsilon_a,epsilon_m,epsilon_s,sigma,omega_a,omega_m,omega_s,num_comps,n_countries);

EV_open_economy = calc_EV(C_vec_open_economy_counterfactual, income_per_capita_baseline, P_a_vec_baseline, P_m_vec_baseline, P_s_vec_baseline, ...
    epsilon_a,epsilon_m,epsilon_s,sigma,omega_a,omega_m,omega_s,num_comps,n_countries);

% Save matrices and output results

GDP_Tornqvist_mat = [GDP_Tornqvist_baseline' GDP_Tornqvist_no_reallocation' GDP_Tornqvist_closed_economy' GDP_Tornqvist_open_economy'];
GDP_Tornqvist_changes_mat = [GDP_Tornqvist_baseline'./GDP_Tornqvist_baseline' ...
    GDP_Tornqvist_no_reallocation'./GDP_Tornqvist_baseline' ...
    GDP_Tornqvist_closed_economy'./GDP_Tornqvist_baseline' ...
    GDP_Tornqvist_open_economy'./GDP_Tornqvist_baseline'];

ag_GDP_share_mat = [ag_GDP_share_baseline' ag_GDP_share_closed_economy_counterfactual' ag_GDP_share_open_economy_counterfactual'];
manu_GDP_share_mat = [manu_GDP_share_baseline' manu_GDP_share_closed_economy_counterfactual' manu_GDP_share_open_economy_counterfactual'];
serv_GDP_share_mat = [serv_GDP_share_baseline' serv_GDP_share_closed_economy_counterfactual' serv_GDP_share_open_economy_counterfactual'];

ag_domestic_expenditures_mat = [ag_domestic_expenditures_vec_baseline ag_domestic_expenditures_vec_open_economy];
manu_domestic_expenditures_mat = [manu_domestic_expenditures_vec_baseline manu_domestic_expenditures_vec_open_economy];

EV_mat = [EV_naive' EV_closed_economy' EV_open_economy'];

food_prices_mat = [P_a_vec_baseline' P_a_naive_counterfactual_vec' P_a_vec_open_economy_counterfactual'];

Rawlsian_losses = Rawlsian_CE(L_vec,P_a_vec_baseline,P_m_vec_baseline,P_s_vec_baseline,...
    C_vec_baseline,C_vec_open_economy_counterfactual,...
    epsilon_a,epsilon_m,epsilon_s,sigma,omega_a,omega_m,omega_s,income_per_capita_baseline);

Rawlsian_losses_log = Rawlsian_CE_log(L_vec,income_per_capita_baseline,EV_open_economy);

Rawlsian_losses_CRRA = Rawlsian_CE_CRRA(L_vec,income_per_capita_baseline,EV_open_economy);

outputs = [GDP_Tornqvist_changes_mat EV_mat ag_GDP_share_mat manu_GDP_share_mat serv_GDP_share_mat ...
    ag_domestic_expenditures_mat ag_productivity_shocks' manu_shocks' serv_shocks' ...
    ag_net_exports_baseline' ag_net_exports_open_economy_counterfactual' food_prices_mat ...
    Rawlsian_losses_log.*ones(n_countries,1) Rawlsian_losses_CRRA.*ones(n_countries,1) Rawlsian_losses.*ones(n_countries,1) ...
    manu_net_exports_baseline' manu_net_exports_open_economy_counterfactual' ...
    manu_domestic_expenditures_mat];

%dlmwrite('results_main_counterfactuals.csv',outputs)



